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We present computer simulations of a simple bead-spring model for polymer melts with in- 
tramolecular barriers. By systematically tuning the strength of the barriers, we investigate their 
role on the glass transition. Dynamic observables are analyzed within the framework of the Mode 
Coupling Theory (MCT). Critical nonergodicity parameters, critical temperatures and dynamic 
exponents are obtained from consistent fits of simulation data to MCT asymptotic laws. The so- 
obtained MCT A-exponent increases from standard values for fully-flexible chains to values close to 
the upper limit for stiff chains. In analogy with systems exhibiting higher-order MCT transitions, we 
suggest that the observed large A-values arise form the interplay between two distinct mechanisms 
for dynamic arrest: general packing effects and polymer-specific intramolecular barriers. We com- 
pare simulation results with numerical solutions of the MCT equations for polymer systems, within 
the polymer reference interaction site model (PRISM) for static correlations. We verify that the 
approximations introduced by the PRISM are fulfilled by simulations, with the same quality for all 
the range of investigated barrier strength. The numerical solutions reproduce the qualitative trends 
of simulations for the dependence of the nonergodicity parameters and critical temperatures on the 
barrier strength. In particular, the increase of the barrier strength at fixed density increases the 
localization length and the critical temperature. However the qualitative agreement between theory 
and simulation breaks in the limit of stiff chains. We discuss the possible origin of this feature. 

PACS numbers: 64.70.pj, 64.70.qj, 61.20.Ja 



I. INTRODUCTION 

Since they do not easily crystallize, polymers are prob- 
ably the most extensively studied systems in relation with 
the glass transition phenomenon. Having said this, their 
macromolecular character, and in particular chain con- 
nectivity, must not be forgotten. The most evident ef- 
fect of chain connectivity is the sublinear increase of the 
mean squared displacement (Rouse-like) [l| arising af- 
ter the decaging process, in contrast to the linear regime 
found in non-polymeric glass-formers. Moreover, in the 
case of strongly entangled polymer chains, the reptation 
model predicts other two sublinear regimes between the 
Rouse and linear regimes P, 0, 0] . 

Another particular ingredient of polymers is that, 
apart from fast librations or methyl group rotations 
every motion involves jumps over carbon-carbon rota- 
tional barriers and/or chain conformational changes. In- 
tramolecular barriers play a decisive role in the physical 
properties of polymer systems. Thus, they are respon- 
sible of partial or total crystallization 0, Q . They also 
enhance dynamic features which are usually associated 
to reptation 0, [1] , which controls rheological properties 
Q. Models for semiflexible polymers are of great interest, 
since they can be applied to many important biopolymers 
such as proteins, DNA, rodlike viruses, or actin filaments 
HO) El- Moreover, chain stiffness seems to play an 



important role in the absorption behavior of polymers at 
interfaces [l^, [13 ■ Thus, an understandig of the role of 
intramolecular barriers on structural, dynamic and rheo- 
logical properties of polymers is of practical as well as of 
fundamental interest. 

In this work we investigate, by means of molecular dy- 
namics simulations, the role of intramolecular barriers on 
the glass transition of polymer melts, by systematically 
tuning barrier strength in a simple bead-spring model. 
We discuss the obtained results within the framework 
of the Mode Coupling Theory (MCT) of the glass tran- 
sition [13, [H, [II El [S We extend preliminary 
results reported by us in Ref. Q by testing a large set 
of predictions, including the factorization theorem and 
time-temperature superposition principle. A consistent 
set of dynamic exponents associated to asymptotic scal- 
ing laws is obtained. By increasing the barrier strength 
a crossover is observed for the values of the so-called A- 
exponent. In the limit of fully- flexible chains A takes val- 
ues ~ 0.7, characteristic of simple fluids dominated by 
packing effects. On the contrary, for strong intramolecu- 
lar barriers the A-values approach the upper limit A = 1 
characteristic of higher-order MCT transitions. The lat- 
ter arise in systems with different coexisting mechanisms 
for dynamic arrest [111, [l^, HI] . In the system investigated 
here, the obtained results suggest an interplay between 
general packing effects and polymer-specific intramolec- 
ular barriers. 
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Chong and co-workers [2J, |25| have recently presented 
an extension of the MCT to simple fully-flexible bead- 
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spring models of polymer systems, in the framework of 
the polymer reference interaction site model (PRISM) 
p^ . I27I . [2^ . In this formalism each molecule is divided 
into interaction sites corresponding to monomers. A key 
assumption of the PRISM is the replacement of the site- 
specific intermolecular surroundings of a monomer by 
an averaged one (equivalent site approximation), while 
keeping the fully intramolecular dependence. We have 
tested the PRISM approximations used by MCT in the 
polymer model here investigated, which incorporates in- 
tramolecular barriers. Likewise, we have solved the MCT 
equations for the location of the MCT 'glass transition' 
temperatures (MCT critical temperatures) and for the 
nonergodicity parameters, which quantify the stability 
of density fluctuations in the reciprocal space. We com- 
pare solutions of the MCT equations with the results 
obtained from the phenomenological analysis of the sim- 
ulation data. We observe that the theory reproduces 
qualitative trends in the nonergodicity parameters and 
critical temperatures. However, the agreement breaks as 
the limit of stiff chains is approached. We discuss the 
possible origins of this feature. 

The article is organized as follows. In Section II we 
describe the model and give simulation details. Static 
correlators are shown in Section III. Moreover the PRISM 
approximations are tested for representative values of the 
barrier strength. Section IV presents qualitative dynamic 
trends as a function of the barrier strength. In Section V 
we summarize the universal predictions of the MCT and 
the equations of motion of the version for polymer melts 
introduced by Chong and co-workers. In Section VI we 
perform a phenomenological analysis of simulation data 
within the MCT, by testing universal scaling laws and 
deriving their associated dynamic exponents. In Section 
VII we compare the results of the former analysis with 
numerical solutions of the MCT equations. We discuss 
the observed differences for stiff chains in Section VIII. 
Conclusions are given in Section IX. 

II. MODEL AND SIMULATION DETAILS 

We have performed molecular dynamics (MD) simu- 
lations of a bead-spring model for which we have im- 
plemented bending and torsional intramolecular barri- 
ers. The monomer- monomer interaction is given by a 
corrected soft-sphere potential 

V{r)=At[{a/rf^ -C^^C2{rlafl (1) 

where e = 1 and a = \. The potential V(r) is set to 
zero beyond the cutoff distance r > ca, with c = 1.15. 
The values Co = 7c^^^ and C2 = 6c~^^ guarantee con- 
tinuity of potential and forces at r = ccr. The potential 
V{r) is purely repulsive. It does not show local min- 
ima within the interaction range r < ca. Thus, it drives 
dynamic arrest only through packing effects. Along the 
chain backbone, of N monomers, an additional finitely- 
extensible nonlinear elastic (FENE) potential [2^ [S^ is 
used to introduce bonds between consecutive monomers: 

l^FENE(r) = -eKFRlHl - {Roa)-^r% (2) 



where Kp = 15 and Rq = 1.5. The superposition of 
potentials H]) and ^ provides an effective bond potential 
for consecutive monomers with a sharp minimum at r « 

0. 985, which makes bond crossing impossible. 
Intramolecular barriers are implemented by means of 

a combined bending Vb, and torsional potential Vr- We 
have used the potentials proposed by Bulacu and van der 
Giessen in Refs. [S, JlJ. The bending potential acts on 
three consecutive monomers along the chain. The angle 
between adjacent pairs of bonds is mantained close to 
the equilibrium value 9o — 109.5° by the cosine harmonic 
bending potential 

VB{e^) = (ei^B/2)(cos0, - cos^o)', (3) 

where 6i is the bending angle between consecutive 
monomers i — 1, i and i + 1 (with 2 < i < — 1). 

The torsional potential constrains the dihedral angle 
which is defined for the consecutive monomers i — 

1, i, i + 1 and i + 2 (with 2 < i < A^ — 2), as the angle 
between the two planes defined by the sets — i, i + 
and {i, i + I, i + 2). The form of this potential is 

VT(6'i, ^'i+i, 0i,i+i) — 
3 

eATx sin^ di sin^ 6'i+i ^ a„ cos" (f>iA+i ■ (4) 

n=0 

The values of the coefficients a„ ar e an — 3.00, ai — 
-5.90, 02 = 2.06, and 03 = 10.95 [1, [Slj. The torsional 
potential depends both on the dihedral angle and 
on the bending angles 9i and As noted in Refs. 

[li mi I numerical instabilities arising when two consec- 
utive bonds align are naturally eliminated by choosing 
the torsional potential without the need of imposing 
rigid constraints on the bending angles. 

In the following, temperature T, time i, distance, wave 
vector q, and monomer density p are given respectively 
in units of e/ks (with ks the Boltzmann constant), 
aijnjt)^!'^ (with m the monomer mass), u, a~^, and 
a^^ . We investigate, at fixed monomer density p — 1.0, 
the temperature dependence of the dynamics for different 
values of the bending and torsion strength, (ATbjATt) — 
(0,0), (4,0.1), (8,0.2), (15,0.5), (25,1), (25,4), and (35,4). 
In the following, all the data presented in the figures and 
discussed in the main text will correspond to p = 1.0. 
This value will not be, in general, explicitly mentioned 
there. We have also studied the case (ATbj^t) = (35,4) 
at density p = 0.93. The specific information of this 
case is given in Table U (see below) . We investigate typ- 
ically 8-10 different temperatures for each set of values 

(/^B,ifT). 

We simulate 300 chains, each chain consisting of A^ = 
10 monomers of mass m = 1, placed in a cubic sim- 
ulation box of lenght Lbox — \^A2'lh for p — 1.0, or 
ibox = 14.7756 for p — 0.93, with periodic boundary 
conditions. Equations of motion are integrated by us- 
ing the velocity Verlet scheme [3^ . Computational ex- 
pense is reduced by implementing a linked-cell method 
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[S^ l- We use a time step ranging from 10~^ to 5 x 10~^. 
We take shorter and longer steps for respectively higher 
and lower values of temperatures, bending and torsional 
constants. The system is prepared by placing and grow- 
ing the chains randomly in the simulation box, with a 
constraint avoiding monomer core overlap. The initial 
monomer density is p = 0.375. Equilibration consists of 
a first run where the box is rescaled periodically by a fac- 
tor 0.99 < / < 1 until the target density p is reached, and 
a second isochoric run at that p. Thermalization at the 
target T is achieved by periodic velocity rescaling. After 
reaching equilibrium, energy, pressure, chain radii of gy- 
ration, and end-to-end distances show no drift. Likewise, 
dynamic correlators show no aging effects. Once the sys- 
tem is equilibrated, a microcanonical run is performed for 
production of configurations, from which static and dy- 
namic correlators are computed. Static correlators pre- 
sented here are averaged over typically 300 equispaced 
configurations. Dynamic correlators are averaged over 
typically 40 equispaced time origins. The typical dura- 
tion of a production run is of 40-200 million time steps 
for respectively high and low temperatures. 

III. STATIC PROPERTIES 

a)Orientational correlations 
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Simulation results presented in this work correspond 
to isotropic phases. We do not observe signatures of 
global orientational order induced by chain stiffness for 
the investigated state points. Thus, by measuring the 
quantity P2{Q) — (3(cos^ 9} — l)/2, where O is the an- 
gle between the end-to-end vectors of two chains, and 
averaging it over all pairs of distinct chains, we obtain 
in all cases values |-P2(0)| < 3 x 10~^. This is illus- 
trated in Fig. [U which shows the time evolution of 
P2(0) along a typical simulation window, both for fully- 
flexible chains, {Kb,Kt) = (0,0), and for representative 
stiff chains, {Kb,Kt) = (35,4). 

Local orientational order is also negligible. This is evi- 
denced by computing a similar correlator -P2(0; ?'cm)- In 
this case the average is performed only over pairs of dis- 
tinct chains for which the distance between their respec- 
tive centers-of-mass is less than rem- Fig. [T] displays, for 
the former cases of fully-flexible and stiff chains, data of 
P2(0;?'cm) for several values of rem- Negligible values of 
^2(0; Tern) are obtained for rem > 2.0. Thus, the time av- 
erage over the simulation time window, tgim, provides val- 
ues |(P2(0;?-cm > 2.0))t,i„| < 0.02. By comparing both 
panels we conclude that chain stiffness does not induce 
a significant increase, if any, of local orientational order. 
Weak local orientational order | (^2(0; ?'cm))tsi„ | ^ 0.1 
is observed only for very small interchain distances (see 
data for P2(0;?'cm — 1-4)). Again, the introduction of 
chain stiffness does not induce clear changes in the orien- 
tational order at this length scale. Note that for small rem 
differences in the represented data for different barrier 
strength may even be statistical artifacts, arising from 



FIG. 1: Time evolution of the global and local orientational 
parameter (see text) for fully-flexible (top) and stiff chains 
with {Kb,Kt) ~ (35,4) (bottom), for two selected low tem- 
peratures. 



the small number of neighboring chains within such dis- 
tances and the limitted time of observation (see the am- 
plitude of the fluctuations in data of Fig. [T]) . 

b) Static structure factors and chain form factors 

Now we present results for static structure factors and 
chain form factors, both for fully-flexible chains and for 
a representative case of stiff chains. Let us consider an 
isotropic homogeneous system of volume V containing n 
identical chains of N monomers. The densities of chains 
and monomers are respectively denoted by pc = n/V and 
p — nN/V. Let us denote the location of a monomer 
along its chain by the index 1 < a < TV. The site-site 
static structure factor for monomers of indices a and b is 
defined as: 

Sabiq)^-{pai-ci,0)pb{q,0)). (5) 
n 

Brackets denote ensemble average. The monomer density 
distribution for wave vector q is given by 

n 
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In this expresion r° is the position vector of the ath 
monomer in the jth chain (1 < j < n). The quantity 
Sabiq) can be spUtted into intrachain and interchain a-b 
correlations: 



Sab{q) = ^abiq) + Pchabiq), 



(7) 



or in matrix form, S{q) = w(g) + pch{q). In Eq. ([7]) 
^ab{q) and hab{q) respectively denote the intrachain and 
interchain correlations between monomers of type a and 
h. By averaging over all the possible pairs (a, b) we ob- 
tain the static correlators S{q), uj{q) and h{q), which are 
related through: 



S{q) = uj{q) + ph{q). 



(8) 



In this expression S{q) is the total static structure fac- 
tor, which equivalently can be obtained as S{q) — 
{nNy^{p{~q,0)p{q,0)), where p{q) = J2a=i P<^il) the 
total monomer density distribution. In Eq. ([8]) the chain 
form factor, uj{q), accounts for all the static intrachain 
correlations, while h{q) accounts for all the static inter- 
chain correlations. 

Fig. [2] (top panel) shows simulation results for S{q) 
as a function of temperature for fully-flexible chains, 
(KbjKt) = (0,0). Data for representative stiff chains, 
{Kb,Kt) — (25,1), are shown in the bottom panel. 
In both cases, no signature of crystallization is present. 
Indeed no sharp Bragg peaks are observed. In both 
cases S{q) shows a maximum at Qmax ~ 7.0. Since 
'S'(qmax) comes from the packing in the first shell around 
a monomer, the latter corresponds to a typical distance 
27r/7.0 « 0.90 in the real space between neighboring 
monomers. On cooling, the peak at gmax ~ 7.0 increases 
in intensity, which is a signature of increasing short-range 
order. 

In Fig. [3]we show, for the former values of {Kb, Kt), 
the corresponding results for the form factors ^((j). We 
note that in the case of fully-flexible chains the form fac- 
tor is nearly independent on temperature. The form fac- 
tor for stiff chains exhibits a certain T-dependence, which 
is however rather weak in comparison with that of S{q). 
The T-dependence of ^(q) becomes more clear at low q- 
values. The way the form factor behaves on lowering the 
temperature is directly connected with the values of the 
mean chain end-to-end radius i?oo- Thus, by decreasing 
temperature from T = 2.0 to T = 0.96, the computed 
Rcc increases from 4.8 to 5.5 for the selected stiff chains. 
This leads, for lower T, to a stronger decay in (^{q) at 
low-g. On the other hand, the value i?ee = 3.6 for the 
fully- flexible chains is almost T-independent, leading to 
a negligible T-dependence of uj{q). 

c)Test of the PRISM approximations 

The MOT for polymer melts developed by Chong and 
co-workers [13, [25l i nvokes several approximations of the 
PRISM theory [27|. In this subsection we summarize 
such approximations and test their validity for all the in- 
vestigated range of barrier strength. The site-site direct 
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FIG. 2: Temperature dependence of the static structure factor 
^(g) for fully-flexible chains (top panel) and for chains with 
barrier strength {Kb,Kt) ~ (25, 1) (bottom panel). 
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FIG. 3: Temperature dependence of the form factor Lo{q) for 
fully-flexible chains (lines) and for chains with barrier strength 
{Kb,Kt) = (25,1) (symbols). For clarity, the inset shows 
results in the range of low-g. Different colors correspond to 
different temperatures, following the legends of Fig. [2l 
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correlation function, Cab{Q), is introduced via the general- 
ized Ornstein-Zernike relation for polyatomic molecules, 
or 'reference interaction site model' (RISM) [s^, 

N 

Kb{q) = ^ ^ax{q)Cxy[q)[^yb{q) + Pchyb[q)], (9) 

in which intramolecular contributions are accounted by 
the form factor terms ujab{q)- By inserting ([7|) in Eq. ([9]), 
Cabiq) is related to Suhiq) and tUabiq) as: 

PcCab{q)^cJ-b\q)-S-,\q). (10) 

Here uj~f^ (q) and S^^ (q) arc the elements of, respectively, 
the matrices vf^^{q) and S^^(q), which are defined as the 
inverses of w(q) and S{q). 

In the equivalent- site approximation (which is exact for 
polymer rings) of the PRISM, chain end effects are ne- 
glected and all sites are treated equivalently for inter- 
chain correlations. Thus, Cab is replaced by the average 
over all (a, 6)-pairs: 

Cab{q) « c{q). (11) 

By introducing this approximation in Eq. ^ and aver- 
aging over all (a, 5)-pairs we find h(q) — uj{q)c{q)[uj{q) + 
ph{q)]. By introducing Eq. ([8]) in the latter expression 
we arrive to the scalar equation 

pciq) = l/uj{q) - l/Siq), (12) 

also known as PRISM equation. 

In Figs. [4] and [5] we test the validity of the equivalent- 
site approximation Cab{q) ~ c{q). We calculate Cab{q) 
and c{q) respectively through Eqs. (fTO)) and p2|) . by us- 
ing the quantities uj~^{q), S~^{q), uj{q), and S{q) as com- 
puted from the simulations. Fig. |4] shows results for the 
fully-flexible case. Data for the case {Kb,Kt) = (25, 1) 
are displayed in Fig. \5\ Both data sets correspond to 
the respective lowest investigated temperatures. We use 
a representation analogous to that of Ref. [51]. Thus, 
top and bottom panels in both figures show the com- 
parison of the averaged c{q) with respectively the matrix 
elements Caa{q) and Cab{q)- The data of Fig. |4]are con- 
sistent with results of Ref. [sj for a similar fully-flexible 
bead-spring model. Data in Fig. [5] constitute new results 
for the case of implemented intramolecular barriers. By 
looking at both figures we conclude that the quality of 
the equivalent-site approximation is not altered by the 
introduction of strong intramolecular barriers. Data in 
Fig. [5]display the same trends as in the fully-flexible case. 
Thus, Cabiq) ^ c(g) is an excellent approximation except 
for correlations involving chain end monomers a = 1 (and 
a = iV by symmetry). The latter show deviations from 
c{q) which are moderate around (/max, this g-range being 
the dominating one in the MCT kernel. 

An additional approximation of the PRISM is the ring 
approximation (which is again exact for polymer rings). 
First we define the quantities Sa{q) = J2b=i ^abiq) and 
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FIG. 4: Test of the equivalent site approximation, Eq. 
for fully-flexible chains at T = 0.50. Top and bottom panels 
compare c{q) with respectively matrix elements Caa(?) and 
Casiq)- The insets enhance the region around the wave vector 
9max for the maximum of the static structure factor S{q). 



Sa^iq) = Y^b^i^abil)- By exploiting the fact that 
for a ring polymer Sa{q) is a-independent, i.e., Sa{q) ~ 

^^"'Ef=i^a(g), wefind 

Saiq) « S{q). (13) 

From the definition of S'^i^{q) and Sb{q) the relation 
X^hLi ^abiq)^b{q) = 1 is exact. By introducing the ring 
approximation Sb{q) ~ S{q) the former relation is trans- 
formed into: 

-^,^S{q). (14) 
Sa [q) 

Fig. [S]shows a test of the ring approximation of Eqs. (fT5|) 
(main panels) and ([T^ (insets). This is done both for 
fully-flexible chains (top panel) and for stiff chains with 
(KbjK'x) — (25,1) (bottom panel). The comparison be- 
tween S{q), Sa{q) and 1/S~^{q) as computed from sim- 
ulations is in general excellent, with the same quality for 
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FIG. 5; Test of the equivalent site approximation, Eq. 
for stiff chains with {Kb, Kt) = (25, 1), at T = 0.96. Top and 
bottom panels compare c{q) with respectively matrix elements 
Caa{q) and Casiq)- The insets enhance the region around the 
wave vector gmax for the maximum of the static structure 
factor S{q). 



FIG. 6: Test of the ring approximation, Eqs. (|13I14|I . Top 
panel: fully-flexible case at T = 0.50. Bottom panel: barrier 
strength {Kb,Kt) = (25,1), at T = 0.96. Main panels and 
insets compare S{q) (symbols) with respectively Sa{q) and 
1/Sa^{q), for the sites a = 1,2 and 5 (lines). 



fully-flexible and stiff chains. Only for the end monomers 
a = 1 (and a — N hy symmetry) significant differences 
between S{q) and \/ S^^{q) are observed around the wave 
vector gmax- 

With all these results we conclude that the approxi- 
mations assumed by the PRISM theory and introduced 
in the MCT equations for polymer melts (see below) are 
fulfilled, with the same quality for all the investigated 
range of barrier strength. 

IV. DYNAMIC PROPERTIES 

In this section we show some phenomenological dy- 
namic features induced by the introduction of intramolec- 
ular barriers in our model. Panels in Fig. [7] show the T- 
dependence of the monomer mean squared displacement 
(MSD) for fuUy-fiexible and representative stiff chains 
with {K-QtKt:) — (25, 1). We observe similar features in 



both cases, but also some differences. After the initial 
ballistic regime, a plateau extends over longer times with 
decreasing temperature. This plateau corresponds to the 
caging regime — i.e., the temporary trapping of each 
monomer in the shell of neighboring monomers around it 
— which is usually observed when approaching a liquid- 
glass transition. At longer times, leaving the plateau, a 
crossover to a Rouse-like sublinear regime ((Ar)^) oc t^-^^ 
[sol , im is observed for the fully-flexible case. The fi- 
nal crossover to linear diffusion ((Ar)^) cx t is reached 
at long times only for the highest investigated tempera- 
tures. However, for the case of stiff chains it is difficult 
to discriminate power-law behavior over significant time 
windows. Apparently, the linear diffusive regime is not 
reached within the simulation time window. 

Fig. [8] shows the monomer MSD, for fixed values of 
density p = 1.0 and temperature T ~ 1.5, as a function 
of the barrier strength. Consistently with results in Ref. 
[3l| , we observe that increasing the strenght of the inter- 
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FIG. 7: Temperature dependence of the monomer mean 
squared displacement for fully-flexible (top) and stiff chains 
with {Kb,Kt) = (25, f ) (bottom). The solid and dashed lines 
indicate respectively sublinear (~ t^'^^) and linear behavior. 
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FIG. 9; Temperature dependence of f{q,t) for fully-flexible 
(top) and stiff chains with {Kb,Kt) = (25, 1) (bottom). The 
wave vector is gmax ~ 7, in both cases corresponding to the 
maximum of the static structure factor S{q). 
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FIG. 8: Monomer mean squared displacement, for several 
values of the barrier strength, at fixed density p = 1.0 and 
temperature T = 1.5. 



nal barriers at fixed p and T leads to slower dynamics. 

Fig. [S] shows simulation results at several tempera- 
tures, both for fully-flexible and stiff chains, for the nor- 
malized density-density correlator f{q,t). The latter is 
defined as fiq,t) = {p{-q,0)piq,t))/{p{-q,0)piq,0)}. In 
both cases the correlator is evaluated at the maximum, 



Qmax ~ 7, of the static structure factor S{q). As in the 
case of the MSD, both the fully-flexible and stiff cases 
exhibit the standard behavior in the proximity of a glass 
transition (sol . [s^. After the initial transient regime, 
f{q, t) shows a first decay to a plateau connected with the 
caging regime. On lowering the temperature this plateau 
extends over longer time intervals. At long times, a sec- 
ond decay is observed from the plateau to zero. This 
second decay corresponds to the structural a-relaxation. 

Let us define the relaxation time as a time scale prob- 
ing the a-structural relaxation. This can be done by 
introducing the time for which the correlator for gmax 
takes the value f{quia.x,Tx) = x, provided x is small in 
comparison with the plateau height. Here we use x = 0.2. 
Fig. [TUl shows To. 2 as a funcion of T, for different values 
of the bending and torsional constants. As observed in 
the analysis of the mean squared displacements, increas- 
ing the chain stiffness slows down the dynamics. At fixed 
temperature, the relaxation time for the stiffest investi- 
gated chains increases up to three decades with respect 
to the fuUy-fiexible case. 

In this section we have demonstrated a main dynamic 
feature: the slowing down of the dynamics, at fixed den- 
sity and temperature, by progressively increasing the 
strength of the intramolecular barriers. This feature 
strongly suggests that intramolecular barriers constitute 
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FIG. 10: Temperature dependence of the relaxation times ro.2 
for several values of the barrier strength. 



and additional mechanism for dynamic arrest, coexist- 
ing with the general packing effects induced by density 
and temperature. In the following we summarize the 
main predictions of the Mode Coupling Theory and dis- 
cuss simulation dynamic features within this theoretical 
framework. 

V. MODE COUPLING THEORY: SUMMARY 

In this section we briefly summarize universal dy- 
namic scaling laws concerning the MCT liquid-glass dy- 
namics, and test them in the simulated polymer melt 
for all the investigated range of barrier strength. Ex- 
tensive reviews on MCT can be found, e.g., in Refs. 
[11 M, M, 113, [S E, m, [13. Though initially de- 
rived for simple hard-sphere systems, these predictions 
follow as consequences of the mathematical structure of 
the MCT equations. More specifically, they are associ- 
ated to the bilinear dependence of the memory kernel 
on the density correlators (see below). Thus, MCT pre- 
dicts the same dynamic scaling laws of the monoatomic 
case if such a mathematical structure is retained in sys- 
tems of polyatomic molecules. This is indeed the case of 
the MCT for polymer melts developed by Chong and co- 
workers [2J, [25| (see below). Therefore, the phenomcno- 
logical analysis of our simulation results in terms of MCT 
dynamic scaling laws is justified within the theory. 

By starting from the fundamental Liouville equation of 
motion and using the Mori-Zwanzig projection operator 
formalism one arrives to an integro-differential equation 
for the normalized density-density correlator: 



TnS{q) 



fkBT 
mS{q) 



dt'm{q,t~t')f{q,t') = 0. 



(15) 



This equation is obtained by using projectors 
over the subspace spanned by the densities and 
the longitudinal currents. The memory kernel 



m{q,t — t') (X {Rq{0)Rq{t — t')), where the quanti- 
ties Rg are, within the Mori-Zwanzig formalism, the 
associated fluctuating forces. Since the kernel cannot 
be exactly expressed in terms of f{q,t) and/or its time 
derivatives, Eq. (|15p is not solvable. MCT introduces 
several approximations for the memory kernel, in order 
to provide a closed solvable form of Eq. (fT5|) . These 
approximations are: 

i) It is assumed that the long-time, slow dynamic 
regime of any observable coupled to density fluctuations 
can be expressed as a linear combination of 'mode pairs', 
PkPq-k- Since the exact expression of the correlator of 
the fluctuating forces contains a slow contribution which 
is a linear combination of mode pairs (see e.g., Ref. 
[l7| for details), the former assumption is equivalent 
to neglecting the fast contribution of the fluctuating 
forces. In other words, it is equivalent to assuming a 
large separation between the time scales of the former 
contributions. 

ii) Convolution approximation: three-point static 
correlations are approximated as products of static 
structure factors, 

(p_q(0)pk(0)pq-k(0)) « nNS{q)S{k)S{\q - k|). (16) 

iii) Kawasaki approximation: dynamic four-point corre- 
lations are factorizcd in terms of products of dynamic 
two-point correlations. 



(Pk-q(0)p-q(0)pe_^(t)pe(^)) 



Fiq,t)F{\q-klt), 



(17) 



where the superscript Q denotes evolution with pro- 
jected dynamics (see e.g., Ref. [l3| for details), and 
F{q,t) = {p{~q,0)p{q,t)) are just the unnormalized 
density-density correlators. 

By making use of these three approximations, the 
memory kernel m{q, t) becomes a bilinear form in f{q, t), 



m{q,t) 



(27r)3 



V(q-k,k)/(fc,<)/(|q-k|,t), (18) 



where the vertex V(q — k, k) is given by: 

V(q-k,k) = ^5(g)5(fc)5(|q-k|)x 
[q■kc(k)+q•(q-k)c(|q-k|)]^ 



(19) 



In a monoatomic fluid the direct correlation function c{q) 
is related to the static structure factor via the exact 
Ornstein-Zernike relation [s^l pc{q) = 1 — S~^{q). With 
all this, Eq. ([TS]) has been reduced to a closed set of 
coupled equations which can be solved self-consistently, 
provided S{q) and c{q) are known (the latter are external 
inputs in the MCT equations). 

For the case of systems with molecular architecture, 
Chong and Hirata have obtained [3!|| , by using projectors 
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over site-densities and site-currents, generalized MCT 
equations of motion for site-site correlators Fab{q, t). The 
latter are defined as Fab{q,t) = n^^ {pa{—q,0)pb{q,t)) . 
Note that Fab{q,0) = Sab{q)- The corresponding MCT 
equations of motion are 

Fabiq, t) + Sa'{q)F.b{q, t) 

x=l 

+ V / dt'MaAq,t~t')F,b{q,t') = o, (20) 

where the memory kernel is now given by 



Mabiq,t) = 




{(q • kfcaxik)cbyik)F^y{k, t)Fab{p, t) 
+ (q • k)(q ■ ^)cax{k)cby{p)F^b{k, t)Fay{p,t)], (21) 

with p = q — k. By comparing Eqs. (|20l2ip with Eqs. 
(I15I18I19P we note that the general mathematical struc- 
ture of the kernel (bilinear in site-site correlators), and 
of the MCT equations of motion is retained. 

Except for very small values of N ^ numerical solution of 
the MCT equations for site-site correlators is extremely 
expensive, and further simplifications are needed in or- 
der to obtain a tractable set of equations. For the case 
of simple bead-spring chains, Chong and co-workers have 
reduced [23, '25*1 Eqs. (|20l2ip to a scalar form for f{q,t). 
This is achieved by introducing in Eqs. (|20l2ip the equiv- 
alent site, Eq. pT|) . and ring, Eqs. p3ll4p . approxima- 
tions of the PRISM theory. The so-obtained scalar MCT 
equations of motion, memory kernel and vertex for poly- 
mer chains are formally identical to Eqs. I|15ll8ll9p . The 
polymer character of the system only enters implicitly 
through the PRISM relation pc{q) = l/Lo{q) - l/S{q), 
which differs from the Ornstein-Zernike equation, pc{q) = 
1 — S~^{q), for monoatomic systems. With this, general 
MCT predictions which originate from the mathematical 
structure of Eqs. (|15ll8ll9p will be, due to the mentioned 
formal equivalence, analogous both for monoatomic sys- 
tems and for polymer chains. Now we summarize such 
general predictions. 

In MCT, nonergodic arrested states (glasses) are de- 
fined as those for which density correlators do not exhibit 
full relaxation. More specifically, if we introduce the non- 
ergodicity parameters, defined as fq — limt_>oo /(?, i), 
MCT discriminates between fluid states [fq = 0) and 
glassy states {fq > 0). At the MCT critical tempera- 
ture Tc, the nonergodicity parameters jump from zero to 
nonzero values [40| . In the following we use the notation 
for referring to the critical nonergodicity parameters, 
i.e., the values of fq a.t T — Tc. 

By Laplace transform {t z) of Eqs. (|15ll8p and tak- 
ing the limit z ^ 0, one finds a coupled set of equations 
for the nonergodicity parameters: 

y (2^V(q-k,k)/iq_ki/fc, (22) 



where .^q({/}) denotes a functional, whose explict ex- 
pression is given in the right-hand side of the equation. 
Note that Eq. ([22|) always has the trivial solution fq — 0. 
Thus, glassy states take place when solutions fq>0 also 
exist. 

Given a tagged chain (labeled s), the density distribu- 
tion for the ath monomer of the tagged chain is defined 
as Pa{q) = exp[iq • r^]. The site-site intrachain correla- 
tor is defined as F^^{q, t) = {pl{—q, 0)pl{q, t)). Note that 
F|j(g,0) = uJabiq)- For the derivation of the MCT equa- 
tions for F^^{q, t) we refer to [11]. In this case the reduc- 
tion to a scalar form is not possible. The corresponding 
nonergodicity parameters fabil) = limt^no F^hi q, t) are 
obtained by solving the N x A'^-matrix equation [41| 

N 

rab{<l)= E K.m + ^',]^y^yb{q), (23) 
x,y=l 

with I the identity matrix. The corresponding functional 
^Ibil) is given by 

/(^ 

V^(q-k,k)/,\(fc)/|q^k|, (24) 

with the vertex 
V^(q-k,k) = ^5(|q-k|)[q.(q-k)pc2(|q-k|). (25) 

The normalized self-correlator, usually introduced 
as r(q,t) = ("iV)-^E;=iEli(exp[*q • (r.^t) - 
'"j(O))]): can be equivalently obtained as f^{q,t) = 

^~^J2a=i-FaaiqT^)- Likewise, the corresponding non- 
ergodicity parameters, defined as the long-time limit 
of /"(g,i), can be obtained as /| = T,a=i faai<l)- 
Thus, the solution of Eq. ((23)) also provides trivially the 
nonergodicity parameters for the self-correlator. 

The separation parameter, e = (T — Tc)/Tc, is intro- 
duced to quantify the relative distance to the critical tem- 
perature Tc- We are interested in the behavior of f{q,t) 
in the ergodic fiuid by approaching from above. Thus 
we express the long-time behavior of the density-density 
correlators as: 

f{q,t) = f'q+gq{t), (26) 

where gq{t) quantifies (small) deviations around for 
|e| ^ 0. By introducing Eq. ^ in Eqs. (|15ll8p . ex- 
panding the functional J-'q of Eq. in a power series of 
e|, comparing the so-obtained resulting expressions and 
retaining the lower-order terms (see, e.g., Ref. [s^ for a 
detailed exposition), one finds that gq{t) = hqG{t), where 
hq only depends on q, and G{t) is a q-independent term 
which contains the full time dependence of the deviations 
of fil^t) around Thus, we rewrite Eq. (pS)) as: 

f{q,t)^f^^+hqGit). (27) 
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This expression is known as the first universality of 
the MCT or factorization theorem. It predicts a scahng 
function G{t) (known as the /3-correlator) that is common 
for all the density correlators (since it is g-independent). 
Following the procedure mentioned in the previous para- 
graph [3(5], the function G{t) is found to obey the equa- 
tion: 

a-z^G^{z)^\zL[G^{t)], (28) 

where Giz) and L[G'^{t)] are the Laplace transform of 
respectively G{t) and G^{t). In this equation a = c\e\, 
with c a constant (see [s^ for its explicit expression) , and 
A is another constant given by 



A = ^ e^Ciq, k, |q - k|)efee|q_k| 

qk 



(29) 



The quantities Cq and are respectively the eigenvectors 
of the so-called stability matrix (see below) and its 
traspose, with the normalization conditions e^e^ = 1 
and 6^(1 ~ = 1- The elements of the stability 

matrix are given by 

C^iq,k) = {l-f^,r(^) . (30) 

V Ofk J {^ = /c| 

The terms C'^{q, k, |q — k|) in Eq. ([29]) are given by: 
C^(g,fc,|q-k|) = 



1 / B'^T 



:.(l-/fe)^(l-/|Vk|r 



(31) 



Eq. (|28p for the /3-correlator does not have an ana- 
lytical solution. Still, asymptotic expressions can be ob- 
tained for different time windows. With this idea in mind 
the /3-time scale is first defined as 



(32) 



with io a microscopic time scale and a an exponent. The 
/3-correlator is then rewritten as G{t) = \a\^^^gc{t/Tp). 
By introducing this expression in Eq. ()28|) and taking 
the limits t <^ Tp and t ^ Tp one finds the asymptotic 
solutions IIGII: 



9a{t/Tfi) = [t/Tfi)-^^ 



(33) 



gAt/rp) = -B{t/Tpf t:^Tp, (34) 

where J5 is a constant ^36*1 . The exponents a and b follow 
the constraint 



A = 



r^(i-a) _ r^(i + 6) 

r(l-2a) ~ V{l + 2hy 



(35) 



where T denotes the Euler's Gamma function. According 
to Eqs. (I33I34|) . one finds the asymptotic expressions for 
Eq. (HZl): 

f{q.t)^ f^.+hq{t/to)-^ t^Tp, (36) 



f{q.t) = fl-hq{t/T^f t:>Tf,. (37) 

The analysis of the long-time decay usually includes 
higher-order corrections 36] to Eq. ([57)1 : 

f{q,t) = f^^ - hq{t/T^f + h^^Kt/T^f" + Oit^"). (38) 

The latter is also known as the von Schweidler expansion. 
In these equations Tq is the a-time scale, defined as: 



The exponent 7 follows the constraint: 



1 



1 



7 



2a 2b' 



(39) 



(40) 



Another important prediction of the MCT for states 
approaching from above, is the second universality or 
time-temperature superposition principle (TTSP). This 
prediction arises as a long-time scaling property of the 
MCT equations of motion [ll|. According to the TTSP, 
the long-time decay of any correlator /(q, t) (i.e. the final 
part of the a- relaxation) is invariant under scaling by the 
a-relaxation time Tq,. In other words, for two tempera- 
tures Ti and T2 above one finds 

f{q,t/To.{T,);Ti) = /(g, t/T„(T2); Ta) = /(g,£), (41) 

where f{q,t) is a T-independent master function of the 
normalized time t. While G{t) is common to all correla- 
tors, the master function f{q,t) associated to the TTSP 
is different for each correlator f{q,t). The superposition 
principle implies that the estimated a-relaxation time, 
defined in this work as the time Tx where /((/max, t) takes 
a value x well below the plateau, is proportional to Tq,. 
Thus, it also follows the asymptotic power law 



,(T) cx (T-Te)~ 



(42) 



The a-decay from the plateau to zero is often well 
described by an empirical Kohlrausch- Williams- Watts 
(KWW) function. 



/(g,0=Aexp[-(t/rf)^'] 



(43) 



with Aq, f3q < 1. Note that the latter does not come out 
as an analytical solution of the MCT equations. However 
in the limit q 
that El 



00 of the KWW time t^, MCT predicts 



oc q 



-i/b 



00, 



(44) 



where b is the von Schweidler exponent introduced above. 

The set of equations exposed in this section consti- 
tute a series of universal results which originate from 
the structure of the MCT equations of motion, Eqs. 
(|15ll8ll9p . As mentioned above, the latter were initially 
derived for simple hard-sphere systems, but the corre- 
sponding ones for polymer melts become formally iden- 
tical following the derivation by Chong and co-workers 
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[lliHEl- With this, the scahng laws exposed in this sec- 
tion will also hold in the MCT for polymer melts. Thus, 
the phenomenological analysis of our simulation data in 
terms of such scaling laws is justified within the frame- 
work of MCT. This analysis is presented in the next sec- 
tion. 

VI. MCT ANALYSIS OF SIMULATIONS 

In order to test the factorization theorem, Eq. ^7\ . 
we compute the ratio: 

_ f{q,t)-f{q,t') _ Gjt) - Gjt') 

f{<l,t")-f{q,t') G{t")~G{t') ^^'^ 

where t' and t" are arbitrary times in the /3-regime. 
The ratio for the self-correlators, Rg{t), is defined anal- 
ogously. If the factorization theorem, and then also the 
right-hand side of Eq. (|45p . is fulfilled, the ratios Rq{t) 
and Rg{t) do not depend on the specific correlator. Fig. 
[m shows Rq{t) and Rg{t) over a broad range of wave 
vectors 2.3 < q < 16.5. The data correspond to bar- 
rier strength {Kb,Kt) = (15,0.5) at T = 0.80. The 
fixed times t" = 0.8 and t' = 100 roughly correspond to 
the beginning and the end of the plateau regime. There 
is an intermediate time window of about two decades 
where the data for density-density and self-correlators 
collapse onto a q-independent master curve, while they 
split at both short and late times. Fig. [T2l demonstrates 
that the master curve is, moreover, the same for both 
density-density and self-correlators. Thus, Figs. [11] and 
[T2]demonstrate the validity of the MCT first universality. 

Fig. [13] shows a test of the TTSP, Eq. for the 

density-density correlator evaluated at Qmax (maximun of 
the static structure factor S{q)). The data correspond to 
the case {Kb, Kt) = (15, 0.5) and cover a broad temper- 
ature range 0.80 <T< 1.5. Data collapse onto a master 
curve after rescaling the absolute time by the relaxation 
time To. 2. Thus, the MCT second universality also holds 
for chains with strong intramolecular barriers. 

Solving numerically the MCT equations and determin- 
ing the dynamic exponents (a, 5, 7, A) is in general a dif- 
ficult task. When numerical solutions are not available, 
nonergodicity parameters, prefactors and exponents in 
Eqs. (|36|38|42|44p can be obtained as fit parameters 
from simulation or experimental data (see, e.g., Refs. 
[isl m, [H, |4^, |45|). Consistency of the analysis requires 
that dynamic correlators and relaxation times are de- 
scribed by a common set of exponents, all of them related 
to a single A-parameter through Eqs. ()35|40|) . 

We have performed this consistency test for all the 
investigated range of barrier strength. The following 
figures in this section illustrate, for some representative 
cases, the analysis of simulation data in terms of MCT 
asymptotic laws. Fig. [TJ] shows for a broad g-range 
(2.0 < q < 14.4), fits to the von Schweidler expansion, 
Eq. (|38p (up to second-order terms). Data correspond 
to density-density correlators f{q,t) for the state point 
iKB,KT) = (15,0.5), T = 0.80 (labelled SI), and to 
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FIG. 11: Test of the factorization theorem, Eq. (|45|l . for 
density-density (top) and self-correlators (bottom) at T = 
0.80 and {Kb,Kt) = (15,0.5). The different curves corre- 
spond to equispaced wave vectors in the range 2.3 < q < 16.5. 
The fixed times are t' = 100 and t" = 0.8. 




FIG. 12: Common representation of Rq(t) (symbols) and 
Rq{t) (lines), for selected wavevectors (common colors cor- 
respond to common g- values) . As in Fig. 1111 data correspond 
to r = 0.80 and {Kb,Kt) = (15,0.5), and the selected fixed 
times are t' = 100 and t" = 0.8. 



self-correlators r{q,t) for {Kb,Kt) = (35,4), T = 1.33 
(labelled S2). A good description of the simulation data 
is achieved, for all the range of g-values and over almost 
four time decades, with a fixed ^-exponent (6 — 0.50 and 
0.37 for respectively SI and S2). 

Fig. [15] displays, for the former values of the bar- 
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FIG. 13: Test of the time-temperature superposition princi- 
ple, Eq. pif) . for the density- density correlators at gmax ~ 7, 
for the case {Kb,Kt) = (15,0.5). 




FIG. 14: Symbols: simulations results for density correlators. 
Top panel: f{q,t) for {Kb,Kt) = (15,0.5), at T = 0.80. 
Bottom panel: r{q,t) for {Kb,Kt) = (35,4), at T = 1.33. 
Identical symbols in both panels correspond to identical wave 
vectors q [values are given in panel (a)]. Lines are fits to 
the von Schweidler expansion, Eq. (|3HJ (up to second-order 
terms), with b = 0.50 (top) and 0.37 (bottom). 
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FIG. 15: Critical nonergodicity parameters, as determined 
from fits to Eq. (|38|l . for different barrier strength. Top and 
bottom panels show data for respectively and /|'^. 



rier strength, the g-dependence of the so-obtained criti- 
cal nonergodicity parameters (/^ for f{q,t) and f^'^ for 
f^{q^t)). For comparison, we also include the fuUy- 
fiexible case {Kb,K-t) = (0,0). As deduced from the 
stronger decay of and f^'^ for stronger barriers, the 
introduction of chain stiffness yields a weaker stability 
of density fluctuations. It also induces a weaker lo- 
calization for self-motions at fixed density. Thus, by 
making an approximate fit of f^'^ to Gaussian behavior, 
f^'^ « exp(— we estimate, at fixed p = 1.0, a lo- 
calization length Ic = 0.19, 0.21, and 0.23 for respectively 
{Kb,Kt) = (0,0), (15,0.5), and (35,4). 

Data of self-correlators from the plateau to the limit of 
the simulation window have been fitted to KWW func- 
tions, Eq. (not shown). Fig. [12] shows the q- 
dependence of the so-obtained KWW relaxation times 
for the former values of the barrier strength, at their re- 
spective lowest investigated temperatures. The lines rep- 
resent tests of the MCT prediction oc q~^^'' for large 
q. A good description of the data is obtained with the 
same 6-exponents used for the independently obtained 
von Schweidler fits of Fig. [T31 

Fig. [T71 shows, for the same values of {Kb, Kt) in Fig. 
[TBI a test of the power law ro.2 oc (T — Tc)^'' for the tem- 
perature dependence of the estimated a-relaxation times. 
The fit covers about three time decades. By representing 
the data in terms of the separation parameter T/Tc — 1, 
clearly different 7-exponents are evidenced for different 
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FIG. 16: Symbols: (/-dependence of KWW relaxation times 
for different barrier strength. Lines are fits to cx q~^^*' (see 
text). From top to bottom h = 0.54, 0.50 and 0.37. 
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FIG. 17: Symbols: T-dependence of relaxation times T0.2 for 
different barrier strength. Lines are fits to oc (T — Tc)~^ (see 
text). From top to bottom 7 = 2.60, 2.74 and 3.43. 



barrier strength. A good description of the data is ob- 
tained with the 7- values derived, through Eqs. p5l l40)l . 
from the &- values used in Figs. [M] and [TB] This resuh 
demonstrates the consistency of the MCT analysis for 
the representative examples showed here, which cover all 
the range of investigated barrier strength between fuUy- 
flexible and stiff chains. 

Similar consistent tests (not shown) have been per- 
formed for the rest of investigated systems. Table [J dis- 
plays the results for the so-obtained A-exponents and 
critical temperatures Tc as a function of {K-b.Kt). We 
also include the corresponding value of the mean end- 
to-end radius (computed at T^), which provides a qual- 
itative characterization of chain stiffness. From the nu- 
merical values in Table U a clear correlation between the 
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35 


4 


6.5 


1.23 


0.862 


0.93 


35 


4 


6.9 


1.02 


0.885 



TABLE I: Values of the MCT A-exponents and critical tem- 
peratures Tc for different p and barrier strength. Also in- 
cluded are the mean chain end-to-end radii i?L at Tc- 



strength of the internal barriers and the values of Tc 
and A is unambiguously demonstrated. The interplay be- 
tween packing effects and intramolecular barriers induces 
a progressive increase of Tc at fixed density. A similar ef- 
fect is observed for the A-exponent, which increases from 
A = 0.761 for fully-flexible chains to A = 0.885 for the 
stiffest investigated chains. The smallest A-values in Ta- 
ble U are typical of simple glass-formers as the archetype 
hard-sphere fluid (A = 0.74) (36j where dynamic arrest is 
driven by packing effects. The largest ones, A < 0.9, are 
similar to those observed in realistic models of polymer 
melts which incorporate the chemical structure of the 
chains. Some examples include poly(vinyl mcthylcther) 
(4^ . polybutadiene [i^l or poly (vinyl ethylene) ^48], with 
respective values of A = 0.87, 0.93, and 0.93. 

Thus, the analysis presented here rationalizes the dif- 
ference in the MCT exponents between fully-flexible 
bead-spring models and real polymers. The systematic 
study performed by tuning the barrier strength suggests 
that large A-exponents in real polymers arise from the 
interplay between two distinct mechanisms for dynamic 
arrest. These are general packing effects and polymer- 
specific intramolecular barriers. Large A-values arising 
from the interplay between distinct arrest mechanisms 
have been observed in systems of very different nature, 
as short-ranged attractive colloids [13, IH, US] (competi- 
tion between hard-sphere repulsion and short-ranged re- 
versible bonding) , polymer blends [5l|, and colloidal 
mixtures with strong dynamic asymmetry [ssl . [s^ l (bulk- 
like caging and matrix-induced confinement), or densified 
silica [5 5] (presumably bonding and packing). Numerical 
solutions of the MCT equations in short-ranged attrac- 
tive colloids and quenched-annealed mixtures jl^ 
have revealed the existence of higher-order MCT transi- 
tions, which are characterized by the upper limit A = 1. 
Whether higher-order MCT transitions are present at 
some region of the control parameter space of the in- 
vestigated model is an open question. 

In this section we have performed a phenomenological 
analysis of the simulation data within the framework of 
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MCT. In the next section the observed trends are com- 
pared with numerical solutions of the MCT equations. 

VII. SOLUTION OF THE MCT EQUATIONS 



We have solved Eqs. (|^^ and for the nonergod- 
icity parameters, for all the investigated range of barrier 
strength. In analogy with the procedure exposed in, e.g., 
Refs. [H, 1131, the integrals over the reciprocal space in 
the corresponding MCT functionals of Eqs. ([^ and (IM)) 
are discretised to a grid of M = 600 equispaced points, 
with g-spacing l^q — 0.1, leading to the expressions: 



1-/. 



X [{xl + ; 



xl)c{k) + (x\ +xl- xl)c{p)fhU (46) 



and 



E E ^-^S{p)[{xl +xl^ xl)c{p)rmk)U. (47) 



k -^p 



In these expressions the wavevectors are defined as q — 



XfcAg, and p 



XpAq, 



with 



XqAq, k 

1/2, 3/2, ...1199/2. The prime at the sums over Xp means 
that the latter are restricted to Xp-values following the 
condition \xq — Xk\ + l/'2 < Xp < Xq + Xk ~ 1/2. 

The solutions of Eq. PS|) are found by a standard it- 
erative procedure /^'''^/[l — /|^^] = ^qilf''}), with j 
the iteration step, and with the initial condition — 1. 
It can be demonstrated that the stability matrix in Eq. 
pop has always a maximum non-degenerate eigenvalue 
E < 1, which takes the upper value _Ec = 1 at the crit- 
ical point [3^. Thus, by following the drift of E with 
changing temperature it is possible to bracket the values 
of the critical nonergodicity parameters and the crit- 
ical temperature Tc, with very high precision. Once the 
values of are obtained, they are fixed in the functional 
of Eq. (|47| . and a small number of iterations is needed 
to find the corresponding critical values falil)- Finally, 
the critical nonergodicity parameters for self-correlations 
are obtained as f^'' = N-^ ^f^^ falil)- 

Following the procedure exposed above, we solved Eq. 
(US)) by inserting as external inputs the structural quan- 
tities, S{q) and c{q), as directly computed from the sim- 
ulations. However, as previously reported in Ref. (25| 
for fully-flexible chains, a MCT transition was not ob- 
served for any of the investigated barrier strength. This 
means that the theoretical critical temperature Tc is be- 
low the lowest simulation temperature for which equili- 
bration was possible. This result is different from the 
usual observation in non-polymeric systems, for which 
the theoretical critical point is accessible in simulation 
time scales. The reason of this difference is, in some 
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FIG. 18: Comparison of the critical nonergodicity parameters 
fq (top) and fq'^ (bottom) obtained from MCT calculations 
(solid lines) with those obtained from the phenomenological 
analysis of the simulation data (symbols, dashed lines are 
guides for the eyes). Results are shown both for fully-flexible 
chains and for representative stifi' chains with {Kb,Kt) = 
(15,0.5). 



way, related with the unability to crystallize of bead- 
spring models, which avoids a fast growing of peaks un- 
der cooling in the static structure factor S{q), leading to 
MCT kernels which are not sufficiently strong to provide 
nonzero solutions of fq. 

Since static correlations computed from our equilib- 
rium simulations do not induce a MCT transition, we are 
forced to use a structural theory for estimating S{q) and 
c(g) at lower temperatures, which will allows us to insert 
them in the MCT equations and to search for the criti- 
cal temperature. Thus, we solve numerically the PRISM 
equation 



pc{q) = l/uj{q) - 1/S{q), 



(48) 



with the Percus-Yevick (PY) closure relation [3^ for the 
non-bonded potential V{r) of Eq. ([T|). The PY relation 
is given by: 

c(r) = [1 - exp(y(r)/fcsT)](/i(r) + 1), (49) 
where c(r) and h(r) are the Fourier transforms in the 



15 



1.2 
T 

c 

0.9 




0.6 



0.3, 



"1 ' I •~ 

K,,=35 K,,,=4 
" ' O 



Kg=25 K.|,=l 
O 



Kg=I5 K.|,=().5 
O 



Kg=4 Kj.=0.1 ^ 

° n O 

■Kb=o V o ° 
o 

□ o 



7 

ee 



FIG. 19; Critical temperature Tc as a function of the end- 
to-end radius R^^. The theoretical T^'^'^^, as obtained by 
solving the MCT equations, is compared with the from 
the phenomenological analysis of the simulation data. 



Kt < 0.5, the dependence of T^^'^ on R'^^ roughly dis- 
plays the same slope as for T^^, with a shift factor 
yMD/yMCT ^ 2 25. Similar shifts between simulation 
and theory, which have their origin in the mean-field 
character of the MCT, are observed in other systems 
pS , [stI . Issj . The range of barrier strength for which 
jiMCT g^j^j^ T^f^^ are roughly parallel is significant. Note 
that for {Kb, Kt) = (8, 0.2) the end-to-end radius R^^ is 
a factor 1.3 longer than for fully-flexible chains. 

By further increasing chain stiffness the differences be- 
tween and T^'~^'^ progressively increase. We ob- 
serve a saturation of the theoretical T^'~^'^ around « 0.55, 
while the simulation T^^ grows up to a value of 1.23 
for the stiffest investigated chains. Thus, the agreement 
between theory and simulation clearly breaks for stiff 
chains. 

Finally, we have computed the corresponding theo- 
retical A-exponents according to the definitions of Eqs. 
(|29l30l3ip . We find an almost constant value of A « 0.72 
for all the investigated range of barrier strength. This 
result is clearly different from the observations in the 
phenomenological MCT analysis of simulation data (see 
data in Table [l| , which provides a strong dependence of 
A on the barrier strength. 

VIII. DISCUSSION 



real space of c{q) and h{q). The coupled set of nonlinear 
equations (I48I49P is solved by a standard Picard iteration 
method [56] for the quantity T{r) = h{r) — c{r), which is 
a smooth function over all the range of r. The form factor 
w{q) is an external input in this procedure. We observed 
(see above) that uj{q) exhibits a very weak temperature 
dependence in comparison to the total static structure 
factor S{q). Thus we just use for each barrier strength 
the i^iq), as computed from the simulations, at the lowest 
temperature for which equilibration was possible. 

In Fig. [18] we show a comparison of the critical non- 
ergodicity parameters (top panel) and /^'^ (bottom 
panel) as obtained from numerical solution of the MCT 
equations, with the results of the fitting procedure of sim- 
ulation data (see above). The theoretical results qualita- 
tively reproduce the simulation trends, and in particular 
the observation that at fixed density the intramolecular 
barriers induce a weaker localization length. Quantita- 
tively, the MCT solutions overstimate the amplitude of 
the nonergodicity parameters, except in the low-q region 
of for which MCT clearly understimates the results. 

In figure we show a representation of the critical 
temperature Tc as a function of the end-to-end radius 
R^^, which quantifies chain stiffness. Values of Tc ob- 
tained from the phenomenological analysis of the simu- 
lations (Tc'^'^) and from the numerical solutions of the 
MCT equations (Tc'^^'^) are compared. We note that 
Tc'^^ seems to grow monotonously with chain stiffness. 
This trend is well reproduced by the theory for low and 
moderate values of the internal barriers. Thus, for val- 
ues of bending and torsional constants Kb < 15 and 



The results reported in the previous section show that, 
though reproducing some qualitative simulation trends 
for low and moderate barriers, numerical solutions of the 
MCT equations exhibit important differences with sim- 
ulation values as the limit of stiff chains is approached. 
Another result to be understood is the clear disagree- 
ment between the almost constant value of the theoret- 
ical A-exponent, and the observed strong dependence of 
simulation values on the barrier strength. 

The observed disagreement between theory and simu- 
lation for strong barriers does not seem to be related with 
the failure of the PRISM approximations for stiff chains, 
which have been introduced in the derivation of the MCT 
equations. Indeed we have shown that the quality of the 
used PRISM approximations is the same for fully-flexible 
and stiff chains (Figs. IH [5]and[n]). Having said this, it 
might be argued that the theory is simply wrong: the 
phenomenological MCT analysis is apparently success- 
ful, but one finds that it has little to do with the theory, 
for stiff chains, when solving the MCT equations. How- 
ever, we remind that the phenomenological analysis has 
shown, for all the investigated range of barrier strength: 
i) the validity of the two MCT universalities, i.e., the 
factorization theorem (Figs. \TI\ and [H]) and the TTSP 
(Fig. [T3| ii) the possibility of a good description of dif- 
ferent dynamic observables (Figs. [HI [TBI and [T7)) with a 
set of dynamic exponents which are consistently trans- 
formed, through Eqs. (|35l40p . to a single A-exponent. 

We believe that all these observations, for all the inves- 
tigated cases, are not fortuitous. At this point it must be 
noted that the predictions referred to in points i) and ii) 
arise, within MCT, as a consequence of the mathematical 
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structure of the equations of motion, more precisely they 
originate from the biUnear form of the memory kernel. 
The specific values of the numerical solutions clearly de- 
pend on the coefficients of the bilinear products (which 
enter through the vertices of the kernel), but the factor- 
ization theorem, the TTSP, and the asymptotic scaling 
laws are universal properties provided the kernel is bilin- 
ear. Thus, the results of the phenomenological analysis 
suggest that the underlying physics may be connected to 
a bilinear memory kernel, though for high barriers the 
actual coefficients strongly differ from those introduced 
by MCT through the vertices, thus leading to theoretical 
results which strongly differ from simulations. 

In other words, the present results suggest that there 
may be relevant static contributions for the case of stiff 
chains which are missing in the MCT vertices. Thus, the 
inclusion of such contributions will increase the strength 
of the kernel and will induce the theoretical transition at 
higher values of Tc j which might improve the comparison 
between T^'^'^ and T^^ of Fig. [H Recalling the three 
main approximations of MCT, we suggest that the con- 
volution approximation, Eq. (|16p . might break for stiff 
chains. Though possibly it is not the case for intermolec- 
ular contributions, its breakdown for intramolecular con- 
tributions in stiff chains is plausible. It is known that the 
convolution approximation fails when static correlations 
show a strong directionality at near-neighbor distances, 
as for e.g., network-forming liquids as silica [5^. This di- 
rectionality is clearly enhanced for intrachain correlations 
by increasing the barrier strength, as evidenced by the 
progressively larger values of the end-to-end radius (see 
Table ^ . For the case of silica, it has been shown that 
the explicit inclusion of three-point static correlations in 
the MCT vertex improves significantly the quality of the 
comparison between theory and simulations . A simi- 
lar improvement might be achieved in the present case by 
similarly incorporating the intrachain three-point static 
contributions. Work in this direction is in progress. 

IX. CONCLUSIONS 

We have performed simulations on a simple bead- 
spring model for polymer melts with intramolecular bar- 
riers. The role of such barriers on the glass transition 
has been investigated by systematically tuning the bar- 
rier strength. Dynamic correlators probing the structural 
relaxation have been analyzed in the framework of the 
Mode Coupling Theory. We have obtained critical noner- 
godicity parameters, critical temperatures and dynamic 
exponents of the theory from consistent fits of simula- 



tion data to MCT asymptotic laws. From the analysis 
of the critical nonergodicity parameters we deduce that 
the presence of the barriers induces a weaker localization 
length in the system at fixed density. The increase of the 
barrier strength at fixed density also induces a higher 
critical temperature Tc- The values of the dynamic ex- 
ponents, as obtained from the phenomenological analy- 
sis of the simulation data, exhibit significant differences 
between the limit of fully-flexible and stiff chains. In 
particular the so-called A-exponent takes standard values 
A ~ 0.7 for the fully- flexible case and values approaching 
the upper limit A = 1 for strong intramolecular barriers. 
While the former A-values are characteristic of simple 
systems dominated by packing effects, transitions with 
A « 1 arise in systems with different competing mecha- 
nisms for dynamic arrest. In our systems these large A- 
values suggest a competition between two distinct mech- 
anisms: general packing effects and polymer-specific in- 
tramolecular barriers. 

For a comparison between simulation and theory, we 
have numerically solved the MCT equations, following a 
recent extension of the MCT by Chong and co-workers 
for polymer melts. The approximations assumed by the 
structural PRISM theory, which are introduced in the 
MCT equations, are fulfilled for all the investigated val- 
ues of the barrier strength. We have compared the crit- 
ical nonergodicity parameters and critical temperatures 
Tc, as obtained by solving the MCT equations, with the 
corresponding values from the phenomenological analy- 
sis of the simulation data. The theoretical calculations 
qualitatively reproduce the trends observed in the sim- 
ulations for low and moderate barriers. However strong 
discrepancies are observed as the limit of high barriers is 
approached. The reason for such a disagreement possi- 
bly lies in the nature of the approximations made in the 
derivation of the MCT equations. In particular, the con- 
volution approximation for three-point static correlations 
might be unadequate for stiff chains. We suggest that 
a reformulation of MCT equations for polymer melts, 
explicitly including intrachain three-point static corre- 
lations, might lead to a better agreement between simu- 
lations and theory. Work in this direction is in progress. 
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